C.................... RSPRIM.FOR ..................................
C------ This routine evaluates the ionization rates for EUV impact
C------ Modified extensively in February 1993 to include nighttime
C------ light sources instead of using subroutine NYTE.
C------ z=altitude; zox,zn2,zo2,&tn are neut dens & temp arrays sza=
C------ solar zenith angle,
      SUBROUTINE 
     >  PRIMPR(IJ,Z,ZOX,ZN2,ZO2,HE,SZA,TN,JMAX,SETCHI,IVLPS,EFAC)
      IMPLICIT REAL(A-H,O-Z)
      INTEGER IVERT
      REAL EFAC    !.. Solar eclipse factor
      DOUBLE PRECISION Z,ZOX,ZN2,ZO2,HE,SZA,TN,SETCHI,COLUM,OTHPR1
     > ,OTHPR2,OTHPR3,CHI,ZZ,TNJ,XN,COLUMN,TAU,FLUX,HEPLS,FBSBN
     > ,DISN,TAUGAM,FLUXG,ALTG,XNSIGF,DSPECT,GL,CLNITE
     > ,VN2D,VN4S,VNNO,VO1D,OP2D,OP2P,N2D,N4S,NNO
      !.. odd N densities from vertical grid
      COMMON/VODDN/VN2D(401),VN4S(401),VNNO(401),VO1D(401),IVERT
      COMMON/MIB/OP2D(401),OP2P(401),N2D(401),N4S(401),NNO(401)
C----- common to hold the EUV and photoelectron production rates 
      COMMON/EUVPRD/EUVION(3,12,401),PEXCIT(3,12,401),PEPION(3,12,401)
      COMMON/SIGS/ZFLUX(37),SIGABS(3,37),ZLAM(37),SIGION(3,37),
     > TPOT(3,10),NNI(3),LAMAX
      COMMON/SOL/UVFAC(59),EUV
      COMMON/NPRD/COLUM(3,401),OTHPR1(6,401),OTHPR2(6,401),OTHPR3(6,401)
      COMMON/MSIS/AP(7),DEC,ETRAN,BLON,F107,F107A,IDAY,SEC,GL(401)
C....... common suggested by Kent Miller
      COMMON/PRICOM/PROB
      DIMENSION COLUMN(3),XN(3),PROB(3,6,37),ZOX(401),XSNPLS(37)
     > ,ZN2(401),ZO2(401),TN(401),HE(401),SZA(401),Z(401),FNITE(37)
     >  ,CLNITE(3)
      COMMON/BIN37/KMAX,TORR_FLX(37)
      DATA LMAX/0/, F107SV/0.0/, IPROBS/0/
      !.. Fluxes for nighttime ion production in the 37 wavelength bins of
      !.. Torr et al GRL 1979. The fluxes are set to reproduce the production
      !.. rates in Strobel et al. PSS, p1027, 1980. Note that most bins are 
      !.. set to zero and that the Strobel production rates are scaled by 
      !.. FNFAC to stabilize the O+ solution below 200 km. Note also that
      !.. the wavelengths in FNITE go from largest (#3=HI) to smallest.
      DATA FNITE/9E5,0.0,9E5,2*0.0,9E6,13*0.0,3E5,8*0.0,3E5,8*0.0/
      DATA FNFAC/1.0/

      !.. Only need to evaluate EUV fluxes on the first call each time step
      IF(IJ.EQ.1) THEN
        IF(ABS((F107-F107SV)/F107).GT.0.005) THEN
          !..... update UV flux factors
          IF(NINT(UVFAC(58)).EQ.-1.OR.NINT(UVFAC(58)).EQ.-3)
     >      CALL FACEUV(F107,F107A,UVFAC,ZFLUX)
          IF(NINT(UVFAC(58)).EQ.-2.OR.NINT(UVFAC(58)).EQ.-4)
     >      CALL FACSEE(F107,F107A,UVFAC,ZFLUX)
          CALL FACSR(UVFAC,F107,F107A)
          F107SV=F107
        ENDIF

        !.. See if fluxes read from file
        IF(KMAX.EQ.37) THEN
          DO K=1,LMAX
            ZFLUX(K)=TORR_FLX(K)
          ENDDO
        ENDIF

        !.. call params to get solar flux data and cross sections
        CALL PARAMS(0,LMAX)
        !..  find probability for formation of each state  ........
        IF(IPROBS.EQ.0) CALL PROBS(0,PROB,ZLAM,LMAX,NNI)
        IPROBS=1
      ENDIF

      !... initialization of production rates. 1.0E-15 stabilizes 
      !... e density evaluation at low altitudes in CMINOR
      DO 10 IS=1,3
      DO 10 IK=1,12
         EUVION(IS,IK,IJ)=1.0E-15
 10   CONTINUE
C
      DISN=0.0
      DO 687 I=1,6
      OTHPR2(I,IJ)=1.0E-15
 687  OTHPR1(I,IJ)=1.0E-15
C
C........ Nighttime He+ production is calculated and stored. Attenuated to
C........ avoid excess production at low altitudes
      OTHPR1(2,IJ)= 8E-11* EXP(-1.0E-11*ZN2(IJ)) *HE(IJ)
      DO 786 I=1,3
 786  COLUM(I,IJ)=1.0E+25
      TNJ=TN(IJ)
      XN(1)=ZOX(IJ)
      XN(2)=ZO2(IJ)
      XN(3)=ZN2(IJ)
      ZZ=Z(IJ)*1.0E+5
      CHI=SZA(IJ)

C....... determine if sun is below the horizon ...
C---- Must now do calculation for night production - Feb 93
      ALTG=(6371.0+Z(IJ))*SIN(3.1416-CHI)-6371.0
C....      IF(CHI.GT.1.57.AND.ALTG.LT.85.) RETURN
      IF(Z(IJ).GT.1500) RETURN
C
C...... get column densities for scattered light at night  &&&&&&&&
      CALL SCOLUM(IJ,0.0D0,ZZ,TNJ,XN,CLNITE)
C
C...... evaluate the neutral column density  &&&&&&&&
      CALL SCOLUM(IJ,CHI,ZZ,TNJ,XN,COLUMN)
C........ Store the column densities for the 2-Stream program
      COLUM(1,IJ)=COLUMN(1)
      COLUM(2,IJ)=COLUMN(2)
      COLUM(3,IJ)=COLUMN(3)
C
C........ O2 dissociation by Schumann-Runge UV.
C........ OTHPR1(3,IJ)= dissociation rate. OTHPR1(5,IJ)= Energy
      CALL  SCHUMN(IJ,Z(IJ),ZO2(IJ),COLUMN,OTHPR1(3,IJ),OTHPR1(5,IJ),
     >   EFAC)
C
C---- Calculate hv + NO ion. freq. from Lyman-a (Brasseur & Solomon)
C---- OTHPR2(2,IJ) is photodissociation of NO in the SR bands. 
C---- A small night production from scattered light is included. FREQLY
C---- varies with solar activity using Richards et al. 1994 page 8981
C---- LY_a=2.5E11 (Lean), sigi(NO)=2.0E-18 (Brasseur & Solomon page 329)
      DATA O2LYXS,O2SRXS,FREQSR /1.0E-20,1.0E-21,5.0E-6/
       FREQLY=5.0E-7*(1+4.0E-3*(0.5*(F107+F107A)-80.0))
       OTHPR2(1,IJ)=FREQLY*(EXP(-O2LYXS*COLUMN(2))
     >    +0.001*EXP(-O2LYXS*CLNITE(2)))
       OTHPR2(2,IJ)=FREQSR*(EXP(-O2SRXS*COLUMN(2))
     >    +0.001*EXP(-O2SRXS*CLNITE(2)))
C
      !..  wavelength loop begins here  ----------
      !..  TAU, TAUN = optical depth for day, night 
      HEPLS=0.0
      DO 6 L=1,LMAX
      TAU=0.
      TAUN=0.0
      DO 4 I=1,3
      TAUN=TAUN+SIGABS(I,L)*CLNITE(I)
4     TAU=TAU+SIGABS(I,L)*COLUMN(I)
      IF(TAU.GT.70.0) TAU=70.0
      IF(TAUN.GT.70.0) TAUN=70.0
C---- evaluate nighttime flux and daytime flux
      FLUXN=FNFAC*(F107/75.0)*FNITE(L)*EXP(-TAUN)
      FLUX=ZFLUX(L)*EFAC*EXP(-TAU) + FLUXN       !.. EFAC is eclipse factor

C.......... he+ production. He+ X-S  = 0.25 N2  X-S. HEPRDN = nite He+
      IF(ZLAM(L).LT.500.) HEPLS=HEPLS+HE(IJ)*0.25*SIGION(3,L)*FLUX
C
C--- hv + N -> N+ + e. ion. freq. Richards et al. JGR 1994 page 8989
      DATA XSNPLS/6*0.0,.211,10.294,11.171,10.961,11.244,11.323,12.098
     > ,13.265,12.423,11.951,11.212,11.798,11.758,11.778,11.772,11.503
     > ,11.016,10.578,9.556,8.15,8.302,7.298,6.413,6.399,5.192,5.725
     > ,4.787,3.778,2.3,.878,.286/
      RN4S=VN4S(IJ)
      IF(IVERT.LE.0) RN4S=N4S(IJ)
      OTHPR2(3,IJ)=OTHPR2(3,IJ)+XSNPLS(L)*1.0E-18*FLUX*RN4S
C
      IF(ZLAM(L).GE.600.0) THEN
C........ calculation of total euv absorption-ionization.....
      FBSBN=FLUX*(SIGABS(3,L)-SIGION(3,L))*XN(3)
C....... Save energy absorbed in the photodissociative process
      OTHPR1(4,IJ)=OTHPR1(4,IJ)+1.24E+4*FBSBN/ZLAM(L)
C........ production on atomic nitrogen by dissociation
      DISN=DISN+FBSBN
C...      IF(J.EQ.1) WRITE(6,95) L,ZLAM(L),TAU,FLUX,FBSBN,DISN,HEPLS
 95   FORMAT(I4,F9.1,1P,22E9.1)
C........ take into account the large n2 absorption of lyman gamma(972.54)
      IF(NINT(ZLAM(L)).EQ.975) THEN
      TAUGAM=370E-18*COLUMN(3)
      IF(TAUGAM.GT.70.0)TAUGAM=70.0
      FLUXG=UVFAC(34) *0.82E+9 *EXP(-TAUGAM)
      DISN=DISN+FLUXG*370E-18*XN(3)
      ENDIF
      ENDIF


C***** species loop begins here *****
      DO 304 I=1,3
      !.. O2 heating in the EUV wavelength range
      IF(I.EQ.2) OTHPR1(6,IJ)=OTHPR1(6,IJ)+
     >   1.24E+4*FLUX*(SIGABS(2,L)-SIGION(2,L))*XN(I)/ZLAM(L)
      XNSIGF=XN(I)*SIGION(I,L)*FLUX
      K1=NNI(I)
C
C... dspect=# ions formed by w-l l by ionization of k state of species i
      DO 302 K=1,K1
      DSPECT=XNSIGF*PROB(I,K,L)
C......... store ion production rates .....
      EUVION(I,K,IJ)=EUVION(I,K,IJ)+DSPECT
C
C......... calculation of ion heating rate......
      EUVION(1,10,IJ)=EUVION(1,10,IJ)+DSPECT*TPOT(I,K)
C
 302  CONTINUE
 304  CONTINUE
 6    CONTINUE
C
      !..---   wavelength loop ends here   -----------
C
C.........Store UV disoc of N2 2 atoms produced for every dissociation
      OTHPR1(1,IJ)=2.0*DISN
C........ Transfer He+ production to storage
      OTHPR1(2,IJ)=OTHPR1(2,IJ)+HEPLS
C
 777  RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE SCOLUM(J,CHI,Z,TNJ,XN,COLUMN)
      IMPLICIT REAL(A-H,O-Z)
      DOUBLE PRECISION ZG,CHI,Z,TNJ,XN,COLUMN,ALTG,GE,GR,RE,RP,
     >  SH,XP,Y,ERFY2,CHAPFN,RG,HG,XG,EM,F,G,A,B,C,D,GL
      DOUBLE PRECISION GLTM
C++++ this routine evaluates the neutral column density
C++++ see smith & smith jgr 1972 p 3592, subr ambs is the msis
C++++ neutral atmosphere model for z<120 km, chi=solar zenith
C++++ angle, re & ge radius and grav con for earth
      DIMENSION XN(3),COLUMN(3),SN(5),M(3),DG(19),T(2)
      COMMON/MSIS/AP(7),DEC,ETRAN,BLON,F107,F107A,IDAY,SEC,GL(401)
      DATA A,B,C,D,F,G/1.0606963,0.55643831,1.0619896,1.7245609
     >  ,0.56498823,0.06651874/
      DATA SN/0.0,0.0,0.0,0.0,0.0/
        DATA    EM   ,   M(1) , M(2) , M(3) ,  RE   , GE
     1    / 1.662E-24 ,   16. ,  32. ,  28. ,6.357E8, 980/
      DATA T,ALTG,ERFY2/0.0,0.0,0.0D0,0.0D0/
      DATA DG/19*0.0/

      DO I=1,3
        SN(I)=0.0
        COLUMN(I)=1.E+25
      ENDDO

      IF(Z.LT.70*1.0E5) RETURN   !.. below lower boundary of FLIP

      IF(CHI.LT.1.5708) GO  TO 2938      !.. is sza>90.0 degrees

      !.. calculate grazing incidence parameters
      ALTG=(6371.0E5+Z)*SIN(3.1416-CHI)-6371.0E5
      IF(ALTG.GE.85*1.0E5) THEN
        ZG=ALTG*1.E-5
        GLTM=GL(J)
        CALL AMBS(J,ZG,GLTM,DG,T,CHIX,SAT,GLON,GLATD) !..  neutral density at grazing incidence
        SN(1)=DG(2)
        SN(3)=DG(3)
        SN(2)=DG(4)
        TNJ=T(2)
      ELSE
        RETURN
      ENDIF
      !.. sn(1)=o , sn(2)=o2 , sn(3)=n2 , tnj=tn,  gr=gravity, rp=distance
      !.. to pt p, sh=scale height, rg=distance to pt g, hg=scale height at g

 2938 CONTINUE
      GR=GE*(RE/(RE+Z))**2
      RP=RE+Z
      DO 10 I=1,3
      SH=(1.38E-16*TNJ)/(EM*M(I)*GR)
      XP=RP/SH
      Y=SQRT(0.5*XP)*ABS(COS(CHI))
      IF(Y.GT.100.0) 
     >  WRITE(6,100) I,Z/1.0E5,CHI*57.3,TNJ,EM,REAL(M(I)),GR,RP
  100 FORMAT('WARNING, Y IN COLUMN(I) > 100',I4,1P,9E10.2)
      IF(Y.GT.8) ERFY2=F/(G+Y)
      IF(Y.LT.8) ERFY2=(A+B*Y)/(C+D*Y+Y*Y)
    4 IF(CHI.GT.1.5708)GO  TO 2
      CHAPFN=SQRT(0.5*3.1416*XP)*ERFY2

      COLUMN(I)=XN(I)*SH*CHAPFN
        GO TO 10
    2 RG=RP*SIN(3.1416-CHI)
      HG=1.38E-16*TNJ/
     1    (EM*M(I)*980.*(6371.E5/(6371.E5+ALTG))**2)
      XG=RG/HG
      COLUMN(I)=SQRT(0.5*3.1416*XG)*HG*(2.0*SN(I)-XN(I)*ERFY2)
10       CONTINUE
      RETURN
      END
C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
C........ this program determines the cross sections, solar fluxes, and
C........ given in m. torr et al g.r.l 1979 p771, table 2 and 3. but
C........ the longer wavelengths come first
      SUBROUTINE PARAMS(ISW,LMAX)
      COMMON/SIGS/ZFLUX(37),SIGABS(3,37),ZLAM(37),SIGION(3,37),
     > TPOT(3,10),NNI(3),LAMAX
      COMMON/SOL/UVFAC(59),EUV
      DIMENSION  X1(111),X2(111),X3(18),ZLX(37),ZFX(37)
C
      !... ionization potentials for o,o2 and n2 see kirby et al note the
      !... o2 2(pi)u and 2(sigma-)u , and n2 f2(sigma)u pots are guesses
      !... the sixth n2 potential is for dissociation
      DATA X3/13.6,16.9,18.6,28.5,40.,0.0,12.1,16.5,18.2,20.,
     > 25.,0.,15.6,16.7,18.8,25.,29.,37./
      !.... wavelength data. average is taken for bands
      DATA ZLX/1025.,1031.91,1025.72,975.,977.02,925.,875.,825.,775.,
     > 789.36,770.41,765.15,725.,703.36,675.,625.,629.73,609.76,575.,
     > 584.33,554.31,525.,475.,465.22,425.,375.,368.07,325.,303.78,
     > 303.31,275.,284.15,256.3,225.,175.,125.,75./

      !.... absorption cross sections -- o first ,o2, then n2
      DATA X1/5*0.0,1.315,4.554,3.498,5.091,3.749,3.89,4,10.736,11.46
     > ,17.245,13.365,13.4,13.4,13.024,13.09,12.59,12.059,12.127,11.93
     > ,11.496,9.687,9.84,8.693,7.7,7.68,6.461,7.08,6.05,5.202,3.732
     > ,1.839,.73,  1.346,1.0,1.63,21.108,18.73,12.817,8.562,16.631
     > ,22.145,26.668,18.91,20.8,28.535,27.44,21.919,26.017,32.06
     > ,28.07,26.61,22.79,26.04,24.606,23.101,21.91,20.31,18.118
     > ,18.32,17.438,16.81,16.8,14.387,15.79,13.37,10.9,7.509,3.806
     > ,1.316,  3*0.0,50.988,2.24,9.68,20.249,16.992,33.578,16.487
     > ,14.18,120.49,24.662,26.54,31.755,23.339,23.37,22.79,22.787
     > ,22.4,24.13,24.501,23.471,23.16,21.675,16.395,16.91,13.857
     > ,11.7,11.67,10.493,10.9,10.21,8.392,4.958,2.261,0.72/
      !... ionization cross sections 
      DATA X2/5*0.0,1.315,4.554,3.498,5.091,3.749,3.89,4,10.736,11.46
     > ,17.245,13.365,13.4,13.4,13.024,13.09,12.59,12.059,12.127,11.93
     > ,11.496,9.687,9.84,8.693,7.7,7.68,6.461,7.08,6.05,5.202,3.732
     > ,1.839,.73,  .259,0.0,1.05,13.94,15.54,9.374,5.494,6.413,10.597
     > ,10.191,8.47,11.72,23.805,23.75,21.306,24.937,31.1,26.39
     > ,26.61,22.79,26.04,24.606,23.101,21.91,20.31,18.118
     > ,18.32,17.438,16.81,16.8,14.387,15.79,13.37,10.9,7.509,3.806
     > ,1.316,  8*0.0,14.274,8.86,8.5,65.8,15.06,25.48,29.235
     > ,23.339,23.37,22.79,22.787
     > ,22.4,24.13,24.501,23.471,23.16,21.675,16.395,16.91,13.857
     > ,11.7,11.67,10.493,10.9,10.21,8.392,4.958,2.261,0.72/

      IF(ISW.EQ.0.AND.NNI(1).EQ.5) RETURN  !.. only need to calculate once
      NNI(1)=5
      NNI(2)=5
      NNI(3)=6
      LMAX=37
      IF(ISW.NE.0) WRITE(17,95)
 95   FORMAT(/5X,'EUV fluxes, Photoabsorption, and Photoionization ',
     >  'Cross sections',
     > /4X,'I',5X,'lam',5X,'flux',4X,'sigaOX',3X,'sigaO2'
     > ,3X,'sigaN2',3X,'sigiOX',3X,'sigiO2',3X,'sigiN2',3X,'UVfac')
C
      DO 20 L=1,LMAX
      ZLAM(L)=ZLX(L)
      !..- setting up ionization potentials
      IF(L.LE.6)THEN
         TPOT(1,L)=X3(L)
         TPOT(2,L)=X3(6+L)
         TPOT(3,L)=X3(12+L)
      ENDIF
      !..- setting up cross sections
      DO 10 IS=1,3
      IN=LMAX*(IS-1)+L
      SIGABS(IS,L)=X1(IN)*1.0E-18
      SIGION(IS,L)=X2(IN)*1.0E-18
      IF(SIGABS(IS,L).LT.SIGION(IS,L)) SIGABS(IS,L)=SIGION(IS,L)
 10   CONTINUE

      IF(ISW.EQ.0) GO TO 20
      WRITE(17,90) L,ZLAM(L),ZFLUX(L),(SIGABS(I,L),I=1,3)
     > ,(SIGION(I,L),I=1,3),UVFAC(LMAX+1-L)
 20   CONTINUE
C
      IF(ISW.EQ.0) RETURN
      WRITE(17,94)
 94   FORMAT(/5X,' Ionization potentials for O, O2, N2'
     > ,/2X,'4S   2D   2P   4P   2P*  -   X2   a+A  b4   B2   dis  -'
     > ,'  X2   A2   B2   C2   F2   2s')
      WRITE(17,91) ((TPOT(I,J),J=1,6),I=1,3)
C
      RETURN
 90   FORMAT(1X,I4,F9.2,1P,22E9.2)
 91   FORMAT(22F5.1)
      END
C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE PROBS(ISW,PROB,ZLAM,LMAX,NNI)
C.... program for finding branching ratios (probabilities for various ion
C.... and molecular states) of o,o2,n2
C.... ---refs--- m. torr et al grl 1979 page 771, kirby et al atomic data
C.... and nuclear tables 1979 23,page 63
      DIMENSION YO(37,5),PROB(3,6,37),ZLAM(37),NNI(3)
C...... coefficients of o ionization cross sections from torr et al
C..... table 2
      DATA YO/.19,.486,.952,1.311,1.539,1.77,1.628,1.92,1.925,2.259
     > ,2.559,2.523,3.073,3.34,3.394,3.421,3.65,3.92,3.62,3.61,3.88,4.25
     > ,5.128,4.89,6.739,4.0,3.89,3.749,5.091,3.498,4.554,1.315,5*0.0
     > ,.206,.529,1.171,1.762,2.138,2.62,2.325,2.842,2.849,3.446,3.936
     > ,3.883,4.896,5.37,5.459,5.427,5.67,6.02,5.91,6.17,6.29,6.159
     > ,11.453,6.57,3.997,12*0.0, .134,.345,.768,1.144,1.363,1.63
     > ,1.488,1.92,1.925,2.173,2.558,2.422,2.986,3.22,3.274,3.211,3.27
     > ,3.15,3.494,3.62,3.23,2.956,0.664,14*0.0,  .062,.163,.348,.508
     > ,.598,.71,.637,.691,.693,.815,.787,.859,.541,24*0.0, .049,.13
     > ,.278,.366,.412,.35,.383,.307,.308,28*0.0/
C
C....... production of o states from torr et al table 2 (yo array)
C....... need to reverse order of yo to correspond with lambda
      DO 10 L=1,LMAX
      LL=LMAX+1-L
      SUM=YO(LL,1)+YO(LL,2)+YO(LL,3)+YO(LL,4)+YO(LL,5)
      DO 20 I=1,5
      PROB(1,I,L)=0.0
 20   IF(SUM.NE.0.0) PROB(1,I,L)=YO(LL,I)/SUM
 10    CONTINUE
C
C....... call separate subroutines for o2 and n2 probabilities
      DO 30 L=1,LMAX
      CALL PROBO2(1,L,ZLAM(L),PROB,NNI(2))
      CALL PROBN2(1,L,ZLAM(L),PROB,NNI(3))
 30   CONTINUE
C
      IF(ISW.EQ.0) RETURN
      WRITE(17,95)
 95   FORMAT(/5X,' Photoionization branching ratios for O, O2, N2'
     > ,/3X,'Lam    4S   2D   2P   4P   2P*   -   X2   a+A  b4   B2 '
     > ,'  dis   -  X2   A2   B2   C2   F2   2s')
      DO 50 L=1,LMAX
 50   WRITE(17,90) ZLAM(L),((PROB(IS,J,L),J=1,6),IS=1,3)
 90   FORMAT(F8.2,22F5.2)
       RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
         SUBROUTINE PROBN2(ISW,L,ZLAM,PROB,JPTS)
C...... the n2 probabilities are taken from kirby et al tables b and c
C...... the yield of n+ is determined first then the remaining portion
C...... of the cross section is distributed amongst the n2+ ion states
C...... (x,a,b). the dissociation yield is divided between the 3 higher
C...... energy states according to wight et al. j.phys. b, 1976
C...... the 2 other states of kirby et al are not included
C
      DIMENSION A(6),B(6),X(14),Y(14,6),PROB(3,6,37)
      DATA IPTS/14/
      DATA X/50.,210.,240.,280.,300.,332.,428.,500.,600.,660.,660.01,
     > 720.,747.,796./
      DATA Y/5*.32,.3,.46,.404,.308,.308,.308,.42,1.,1.,5*.55,
     > .52,.46,.506,2*.589,.692,.58,2*.0,5*.13,.12,.08,
     > .09,.103,.103,4*0.0,3*0.0,.05,.1,.15,.83,1.,6*.0,3*.0,
     > .3,.4,.79,.17,7*.0,3*1.,.65,.5,.06,8*.0/
C
C...... if zlam is too big set equal to x(max)
      YLAM=ZLAM
      !.. Prevent divide by zero
      IF(ZLAM.GT.X(14)) YLAM=X(14)-1
      IF(ZLAM.LT.X(1)) YLAM=X(1)+1

       YIELD=0.0
C...... determine yield of n+, and store in prob array
       CALL YLDISS(1,YLAM,YIELD)
C
       DO 10 I=1,IPTS
C kjh 6/22/92   NOTE:  I realize the following statement is strange
C   looking, but its purpose is to prevent the CRAY compiler from
C   vectorizing this loop.  (Which it does incorrectly).
      if(i.eq.25)write(6,*)' '
      IF(YLAM.GT.X(I).AND.YLAM.LE.X(I+1))  GO TO 20
 10   CONTINUE
 20   SUM=0.0
C...... fit straight line between points
      DO 30 J=1,JPTS
      A(J)=(Y(I+1,J)-Y(I,J))/(X(I+1)-X(I))
      B(J)=Y(I,J)-A(J)*X(I)
 30   CONTINUE
C...... determine probabilities of n2+ states
      DO 40 J=1,JPTS
      IF(J.LE.3) PROB(3,J,L)=(A(J)*YLAM+B(J))*(1-YIELD)
      IF(J.GT.3) PROB(3,J,L)=(A(J)*YLAM+B(J))*YIELD
      SUM=SUM+PROB(3,J,L)
 40   CONTINUE
C
      IF(SUM.EQ.0.0) RETURN
C....... normalise probabilities
      DO 50 J=1,JPTS
 50   PROB(3,J,L)=PROB(3,J,L)/SUM
      RETURN
      END
C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
         SUBROUTINE YLDISS(ISW,ZLAM,YIELD)
C..... determination of dissociative yield of n+, refer to kirby et al
C..... page 66 and table b
      DIMENSION X(9),Y(9)
      REAL FACTOR
      DATA IPTS/9/
      DATA X/50.,210.,240.,302.,387.,477.,496.,509.,2000./
      DATA Y/.36,.36,.346,.202,.033,.041,.024,0.0,0.0/
C
       DO 10 I=1,IPTS
C kjh 6/22/92   NOTE:  I realize the following statement is strange
C   looking, but its purpose is to prevent the CRAY compiler from
C   vectorizing this loop.  (Which it does incorrectly).
      if(i.eq.25)write(6,*)' '
      IF(ZLAM.GE.X(I).AND.ZLAM.LT.X(I+1))  GO TO 20
 10   CONTINUE
 20   IF(ZLAM.GT.387.AND.ZLAM.LT.477) GO TO 40
      !.. linear interpolation
      FACTOR=(ZLAM-X(I))/(X(I+1)-X(I))
      YIELD=Y(I)+FACTOR*(Y(I+1)-Y(I))
      GO TO 30
C...... parabolic interpolation see formula page 66 kirby et al
 40   YIELD=.0329+8.13E-6*(ZLAM-442)**2
 30   RETURN
      END
C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
         SUBROUTINE PROBO2(ISW,L,ZLAM,PROB,JPTS)
C....... o2 branching ratios are taken from kirby et al table d
C....... columns 4 & 9 are combined. columns 5,6,7&8 are combined
      DIMENSION A(5),B(5),X(20),Y(20,5),PROB(3,6,37)
      DATA IPTS,X/20,304.,323.,454.,461.,504.,537.,556.,573.,584.,598.
     >  ,610.,637.,645.,662.,684.,704.,720.,737.,774.,1026./
      DATA Y/.365,.374,.432,.435,.384,.345,.356,.365,.306,.23,.235,.245,
     > .34,.27,.482,.675,.565,.565,1.,1.,.205,.21,.243,.245,.27,.29,
     > .23,.27,.33,.295,.385,.35,.305,.385,.518,.325,.435,.435,2*.0,
     > .125,.124,.12,.12,.126,.13,.225,.216,.21,.375,.305,.37,.33,.345,
     > 6*.0,.055,.167,.11,.105,.194,.234,.189,.149,.155,.103,.075,.036
     > ,.025,7*0.,.25,.125,.095,.95,.026,15*0./
C
C...... if zlam is too big set equal to x(max)
      YLAM=ZLAM
C...... if zlam is outside range of data values set equal to max or min
      IF(ZLAM.GT.X(20)) YLAM=X(20)
      IF(ZLAM.LE.X(1)) YLAM=X(1)+1.E-3
C
       DO 10 I=1,IPTS
C kjh 6/22/92   NOTE:  I realize the following statement is strange
C   looking, but its purpose is to prevent the CRAY compiler from
C   vectorizing this loop.  (Which it does incorrectly).
      if(i.eq.25)write(6,*)' '
      IF(YLAM.GT.X(I).AND.YLAM.LE.X(I+1))  GO TO 20
 10   CONTINUE
 20   SUM=0.0
C
      DO 30 J=1,JPTS
      A(J)=(Y(I+1,J)-Y(I,J))/(X(I+1)-X(I))
      B(J)=Y(I,J)-A(J)*X(I)
      SUM=SUM+A(J)*YLAM+B(J)
 30   CONTINUE
C
      DO 40 J=1,JPTS
      PROB(2,J,L)=(A(J)*YLAM+B(J))/SUM
 40   CONTINUE
C
      RETURN
      END
C::::::::::::::::::::::::: SCHUMN ::::::::::::::::::::::::::::::::::::::::::::::
C... production of o(1d) by schumann-runge bands
C... The fluxes are from the SEE web site data. 
C... P. Richards March 2011
      SUBROUTINE SCHUMN(J,Z,ZO2,COLUMN,SCHUPR,SCHUHT,EFAC)
      IMPLICIT REAL(A-H,O-Z)
      DOUBLE PRECISION Z,ZO2,SCHUPR,SCHUHT,COLUMN,HSRX,FLD
      COMMON/SOL/UVFAC(59),EUV
      DIMENSION COLUMN(3),SRFLUX(8),SRXS(8),SRLAM(8)
      !.. Schumann-Runge fluxes * 1.0E-11
      DATA SRFLUX/3.3074,1.9699,1.0785,.7083,.5008,.2890,.1651,.1269/
      DATA SRXS/.5,1.5,3.4,6,10,13,15,11/
      DATA SRLAM/1725,1675,1625,1575,1525,1475,1425,1375/

      !.. lmax=# of lambdas in sub. primpr: schuht=heating: schupr=o(1d) prod
      LMAX=37

      DO LSR=1,8
        !.. photoabsorption cross section
        SRXSCT=1.0E-18*SRXS(LSR)
        HSRX=SRXSCT*COLUMN(2)
        IF(HSRX.GT.70)HSRX=70
        !.. attentuated solar flux. EFAC is eclipse factor
        FLD=UVFAC(LMAX+LSR)*1.E+11*SRFLUX(LSR)*EXP(-HSRX)*EFAC
        !.. neutral heating SCHUHT and photodissociation rate SCHUPR
        SCHUHT=SCHUHT+1.24E+4*(FLD*SRXSCT)*ZO2/SRLAM(LSR)
        SCHUPR=SCHUPR+FLD*SRXSCT
        !WRITE(6,90) LSR,SRXSCT,FLD,SCHUPR,UVFAC(LMAX+LSR)*SRFLUX(LSR)
      ENDDO

      SCHUPR=ZO2*SCHUPR
 90   FORMAT(2X,I5,1P,9E10.2)
      JTI=1
      RETURN
      END
C:::::::::::::::::::::::::: SUMPRD :::::::::::::::::::::::::
      SUBROUTINE SUMPRD(JMAX)
C----- common to hold the EUV and photoelectron production rates 
      COMMON/EUVPRD/EUVION(3,12,401),PEXCIT(3,12,401),PEPION(3,12,401)
      COMMON/BPRAUR/PAUION(3,12,401),PAUEXC(3,12,401)
      COMMON/BPRTOT/SUMION(3,12,401),SUMEXC(3,12,401)
C
      DO 30 J=1,JMAX
C----    add contributions from 2 highest O+ metastables to 3 lowest
         EUVION(1,7,J) = EUVION(1,1,J) + EUVION(1,4,J)
         EUVION(1,8,J) = EUVION(1,2,J) + EUVION(1,5,J)/1.3
         EUVION(1,9,J) = EUVION(1,3,J) + EUVION(1,5,J)/4.3
C
         PEPION(1,7,J) = PEPION(1,1,J) + PEPION(1,4,J)
         PEPION(1,8,J) = PEPION(1,2,J) + PEPION(1,5,J)/1.3
         PEPION(1,9,J) = PEPION(1,3,J) + PEPION(1,5,J)/4.3
C
         PAUION(1,7,J) = PAUION(1,1,J) + PAUION(1,4,J)
         PAUION(1,8,J) = PAUION(1,2,J) + PAUION(1,5,J)/1.3
         PAUION(1,9,J) = PAUION(1,3,J) + PAUION(1,5,J)/4.3
 
C----- SUM non-diss. states to get total O2+ and N2+ production
         EUVION(2,7,J)= EUVION(2,1,J)+EUVION(2,2,J)+EUVION(2,3,J)
         PEPION(2,7,J)= PEPION(2,1,J)+PEPION(2,2,J)+PEPION(2,3,J)
         PAUION(2,7,J)= PAUION(2,1,J)+PAUION(2,2,J)+PAUION(2,3,J)
C
         EUVION(3,7,J)= EUVION(3,1,J)+EUVION(3,2,J)+EUVION(3,3,J)
         PEPION(3,7,J)= PEPION(3,1,J)+PEPION(3,2,J)+PEPION(3,3,J)
         PAUION(3,7,J)= PAUION(3,1,J)+PAUION(3,2,J)+PAUION(3,3,J)
 30   CONTINUE
C
      DO 50 I=1,3
      DO 50 K=1,12
      DO 50 J=1,JMAX
        SUMION(I,K,J)=EUVION(I,K,J)+PEPION(I,K,J)+PAUION(I,K,J)
        SUMEXC(I,K,J)=PEXCIT(I,K,J)+PAUEXC(I,K,J)
 50   CONTINUE
C
      RETURN
      END
C::::::::::::::::::::::::::::::: ECFUN ::::::::::::::::::::::::::
C.... This routine calculates the solar eclipse function. It returns the factor
C.... for scaling the solar EUV flux.
C.... 2017-09-08: Redid solar eclipse routine (ECFUN) to split the
C.... irradiance into chromosphere responsible for photoinization
C.... and corona responsible for photoelectrons
C.... Reference in file: Eclipse-Area-of-Overlapping-Circles.doc. 
C.... Written by P. Richards in May 2013. Revised January 2018.
 
      SUBROUTINE ECFUN(SEC,  !.. UT in seconds
     >                IDAY,  !.. FLIP begin day (YYYYddd)
     >                FLAT,  !.. FLIP magnetic latitude to determine eclipse hemisphere
     >                ECDT,  !.. Time step during the eclipse
     >              ECHROM,  !.. Proportion of chromospheric emissions at totality
     >              ECORON)  !.. Proportion of coronal emissions at totality
      IMPLICIT NONE
      INTEGER IDAY
      INTEGER EYYYYDDD       !.. Eclipse year and day from file 
      INTEGER FYYYYDDD       !.. FLIP day (YYYYddd)
      INTEGER EFILE          !.. Eclipse file flag
      INTEGER EFID           !.. Eclipse file unit ID 
      REAL EFAC              !.. Proportion of sun not covered at totality
      REAL ECHROM,ECORON     !.. Proportion of chromospheric and coronal emissions at totality
      REAL CHROM_EMISS       !.. Amount of chromospheric emission at totality (0.0-0.10)
      REAL CORON_EMISS       !.. Amount of coronal emission at totality (0.10-0.30)
      REAL ECDT              !.. Time step during the eclipse
      REAL SEC,SX,UTHRS      !.. UT in seconds, hours
      REAL X,D,AC            !.. factors used to calculate ECFAC
      REAL UTMAX,EDUR        !.. the UT time for maximum eclipse and duration
      REAL FLAT              !.. FLIP latitude
      !.. SMSEPN and SMSEPS are the lateral displacements of the sun and moon centers at 
      !.. time of maximum eclipse in the Northern and Southern hemispheres. SMSEP=0.0 for 
      !.. total eclipse, SMSEP > 2.0 for no eclipse. Specify >2.0 in opposite hemisphere
      REAL SMSEP,SMSEPN,SMSEPS   
      DATA EFID/48/,EFILE/-1/,UTMAX/-99/

      EFAC=1.0
      ECORON=1.0
      ECHROM=1.0
      !... see if eclipse file exists and read data
      IF(EFILE.EQ.0) RETURN
      IF(EFILE.EQ.-1) THEN
        CALL TFILE(EFID,EFILE)
        IF(EFILE.EQ.0) RETURN
  10    READ(EFID,*,ERR=10,END=99) 
     >      EYYYYDDD,UTMAX,EDUR,ECDT,SMSEPN,SMSEPS,
     >      CHROM_EMISS,CORON_EMISS
        CLOSE(EFID)
        IF(SMSEPN.GE.0.AND.SMSEPN.LT.2) 
     >   WRITE(6,'(/A,A,I7,A,F6.2,A,F6.2,A/A,F5.2,A,F5.2)')
     >   ' ** There will be an ECLIPSE in the northern',
     >   ' hemisphere on day=',EYYYYDDD,'; maximizing at UT=',UTMAX,
     >   '; Eclipse time step=',ECDT,' secs;',
     >   '    Chromospheric emission fraction at totality=',CHROM_EMISS,
     >   '; Coronal emission fraction at totality=',CORON_EMISS
        IF(SMSEPS.GE.0.AND.SMSEPS.LT.2) 
     >   WRITE(6,'(/A,A,I7,A,F6.2,A,F6.2,A/A,F5.2,A,F5.2)')
     >   ' ** There will be an ECLIPSE in the southern',
     >   ' hemisphere on day=',EYYYYDDD,'; maximizing at UT=',UTMAX,
     >   ' Eclipse time step=',ECDT,' secs;',
     >   '    Chromospheric emission fraction at totality=',CHROM_EMISS,
     >   '; Coronal emission fraction at totality=',CORON_EMISS
      ENDIF      

      !.. Calculate actual day = FLIP start day + UT elapsed
      CALL ACTUAL_DAY(IDAY,SEC,FYYYYDDD,SX) 
      UTHRS=SX/3600.0

      !.. see if FLIP day is the eclipse day
      IF(FYYYYDDD.NE.EYYYYDDD) RETURN
      !.. See if it is the eclipse time
      IF(ABS(UTHRS-UTMAX).GT.EDUR/2.0) RETURN

      !.. Choose hemisphere for the eclipse
      SMSEP=2.0
      IF(SMSEPN.LT.2.AND.FLAT.GE.0) SMSEP=SMSEPN
      IF(SMSEPS.LT.2.AND.FLAT.LT.0) SMSEP=SMSEPS

      !.. Calculate the attenuation factor. Note that in 
      !.. this calculation the sun and moon have unit diameter
        X=4.0*(UTHRS-UTMAX)/EDUR
        D=0.5*SQRT(X**2+SMSEP**2)    !.. 
        IF(D.GE.1) RETURN
        AC=ACOS(D)-D*SQRT(1-D**2)                 !.. overlapping region
        EFAC=(3.142-2*AC)/3.142                   !.. uncovered fraction = zero at totality
        ECORON=CORON_EMISS+(1-CORON_EMISS)*EFAC   !.. Amount of coronal irradiance at totality
        ECHROM=CHROM_EMISS+(1-CHROM_EMISS)*EFAC   !.. Amount of chromosphere irradiance at totality
 99   CONTINUE
      RETURN
      END
C::::::::::::::::::::::::::::::::: PEFACS :::::::::::::::::::::::::::::::
C...... This function produces the eclipse factors for photoelectron
C...... primary production. Coronal emissions are assumed for most energies. 
C...... Chromosphere emissions are assumed for the 304A fluxes.
C...... Programmed by P. Richards January 2018
      SUBROUTINE PEFACS(ECORON,     !.. IN: coronal fraction for eclipse
     >                  ECHROM,     !.. IN: chromosphere fraction for eclipse
     >                  PECFAC)     !.. OUT: Photoelectron eclipse factors
      IMPLICIT NONE
      INTEGER I
      REAL ECORON,ECHROM,PECFAC(901)
      !.. Fill up the EUV scaling factors for primary photoelectrons
      !.. At first, all are assume corona emissions
      DO I=1,901
        PECFAC(I)=ECORON
      ENDDO
      !.. Now scale 304
      PECFAC(23)=ECHROM   !.. 304 emission
      PECFAC(24)=ECHROM   !.. 304 emission
      PECFAC(25)=ECHROM   !.. 304 emission
      PECFAC(26)=ECHROM   !.. 304 emission
      PECFAC(27)=ECHROM   !.. 304 emission
      !.. WRITE(6,'(A,22F10.2)') ' ECF',(PECFAC(I),I=20,30)
      RETURN
      END
